TODO

  1. Fix pred grid and standardization to account for different sizes!

Intro

In this script, I take the collated stomach data set and calculate aggregates (feeding ratio, total weight of prey groups) and predictor variables for diet data, aggregate to get 1 stomach = 1 row per prey type (not prey individual). I also select only the columns I need for model fitting, join environmental covariates and cpue covariates for cod and flounder, and lastly saduria biomass densities.

Load packages & source functions

# Load libraries, install if needed
library(tidyverse)
#> Warning: package 'tidyr' was built under R version 4.0.5
library(readxl)
library(tidylog)
library(RCurl)
library(RColorBrewer)
#> Warning: package 'RColorBrewer' was built under R version 4.0.5
library(patchwork)
library(janitor)
library(forcats)
library(gapminder)
library(viridis)
library(ggridges)
library(raster)
library(icesDatras)
library(ggalluvial)
library(ggrepel)
library(ncdf4)
library(chron)
library(rnaturalearth)
library(rnaturalearthdata)
library(mapplots)
library(geosphere)
library(quantreg)
#> Warning in .recacheSubclasses(def@className, def, env): undefined subclass
#> "numericVector" of class "Mnumeric"; definition not updated
library(brms)
#> Warning: package 'Rcpp' was built under R version 4.0.5
library(sdmTMB)
options(mc.cores = parallel::detectCores()) 

world <- ne_countries(scale = "medium", returnclass = "sf")

# Source code for map plots
# Source code for map plots
source("/Users/maxlindmark/Dropbox/Max work/R/spatial_metabolic_index/R/functions/map_plot.R")

plot_map_fc <- plot_map_fc + theme(legend.position = "bottom")

theme_set(theme_plot())

# Continuous colors
options(ggplot2.continuous.colour = "viridis")

# Discrete colors
scale_colour_discrete <- function(...) {
  scale_colour_brewer(palette = "Set2")
}

scale_fill_discrete <- function(...) {
  scale_fill_brewer(palette = "Set2")
}

Read data

# Load cache
# qwraps2::lazyload_cache_dir(path = "R/01_sdm_mi_cache/html")

pred_grid1 <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/pred_grid_(1_2).csv")
#> Rows: 129346 Columns: 29
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr  (1): ices_rect
#> dbl (28): X, Y, depth, year, oxy, temp, lon, lat, sub_div, density_saduria, ...
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
pred_grid2 <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/pred_grid_(2_2).csv")
#> Rows: 120107 Columns: 29
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr  (1): ices_rect
#> dbl (28): X, Y, depth, year, oxy, temp, lon, lat, sub_div, density_saduria, ...
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

pred_grid <- bind_rows(pred_grid1, pred_grid2)

pred_grid <- pred_grid %>%
  dplyr::select(-density_saduria, -biomass_spr, -biomass_her, -biomass_spr_sd, -biomass_her_sd, 
                -density_cod, -density_fle, -depth_rec, -temp_rec, -oxy_rec, -density_cod_rec, 
                -density_cod_sd, -density_fle_sd, -density_fle_rec, -density_saduria_rec, 
                -density_saduria_sd, -temp_sd, -oxy_sd) %>% 
  drop_na(oxy, temp)
#> drop_na: no rows removed

cpue <- readr::read_csv("data/full_cpue_length_q1_q4.csv") %>%
  janitor::clean_names() %>% 
  rename(X = x,
         Y = y) %>% 
  filter(quarter %in% c(1, 4))
#> Rows: 1759609 Columns: 18
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr  (6): haul.id, IDx, ices_rect, id_haul_stomach, species, haul.id.size
#> dbl (12): density, year, lat, lon, quarter, sub_div, length_cm, depth, temp,...
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> rename: renamed 2 variables (X, Y)
#> 
#> filter: removed 21,168 rows (1%), 1,738,441 rows remaining

# Match pred grid else we cannot predict later
cpue <- cpue %>% filter(year %in% unique(pred_grid$year))
#> filter: removed 43,092 rows (2%), 1,695,349 rows remaining
pred_grid
#> # A tibble: 249,453 × 11
#>        X     Y depth  year   oxy  temp   lon   lat sub_div ices_rect depth_sd
#>    <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>   <dbl> <chr>        <dbl>
#>  1   306  6012  13.2  1993  4.11  9.78  12.0  54.2      24 37G2            21
#>  2   306  6012  13.2  1994  6.27  9.66  12.0  54.2      24 37G2            21
#>  3   306  6012  13.2  1995  5.42  9.63  12.0  54.2      24 37G2            21
#>  4   306  6012  13.2  1996  6.11  9.99  12.0  54.2      24 37G2            21
#>  5   306  6012  13.2  1997  5.78  9.84  12.0  54.2      24 37G2            21
#>  6   306  6012  13.2  1998  6.45  8.56  12.0  54.2      24 37G2            21
#>  7   306  6012  13.2  1999  5.53 10.2   12.0  54.2      24 37G2            21
#>  8   306  6012  13.2  2000  4.37 12.5   12.0  54.2      24 37G2            21
#>  9   306  6012  13.2  2001  5.81  9.91  12.0  54.2      24 37G2            21
#> 10   306  6012  13.2  2002  5.26  9.92  12.0  54.2      24 37G2            21
#> # … with 249,443 more rows

# Now scale variables in the cpue data by the mean and sd in the prediction grid
# pred_grid$log_depth <- log(pred_grid$depth)
# pred_grid$log_depth2 <- pred_grid$log_depth * pred_grid$log_depth
# mean_depth <- mean(pred_grid$log_depth)
# sd_depth <- sd(pred_grid$log_depth)
# mean_depth2 <- mean(pred_grid$log_depth2)
# sd_depth2 <- sd(pred_grid$log_depth2)

# Calculating a metabolic index

# A_0 is the ratio of constant terms in the Arrhenius equations describing the rates of supply and demand
# B is individual body size (g)
# n is the allometric exponent (difference delta - epsilon)
# E0 is the difference in temperature dependences of supply and demand
# T is temperature in degrees Kelvin
# Tref is an arbitrarily chosen reference temperature (15°C)
# kb is Boltzmann’s constant (eV K−1).

# Parameters from Deutsch et al 2015 Science Supp Fig. S2
A_conc = 9.7*10^-15
E_conc <- 0.72
n <- -0.21
Tref <- 273.15 + 10
kb <- 0.000086173324

pred_grid$T_K <- pred_grid$temp + 273.15
pred_grid$oxy2 <- (pred_grid$oxy * 10^3) / 22.391
pred_grid$phi <- A_conc*(2000^n)*(pred_grid$oxy2 / exp(-E_conc / (kb * pred_grid$T_K)))

# Scale variables in data w.r.t. prediction grid
mean_phi<- mean(pred_grid$phi)
sd_phi <- sd(pred_grid$phi)

mean_depth <- mean(pred_grid$depth)
sd_depth <- sd(pred_grid$depth)

mean_oxy <- mean(pred_grid$oxy)
sd_oxy <- sd(pred_grid$oxy)

mean_temp <- mean(pred_grid$temp)
sd_temp <- sd(pred_grid$temp)

# Calculate phi in data
# Metabolic index in the data
B_small_cod <- 0.01*median(filter(cpue, length_cm <= 30 & species == "cod")$length_cm)^3
#> filter: removed 1,415,326 rows (83%), 280,023 rows remaining
B_large_cod <- 0.01*median(filter(cpue, length_cm > 30 & species == "cod")$length_cm)^3
#> filter: removed 819,148 rows (48%), 876,201 rows remaining

# https://github.com/fate-spatialindicators/SDM_O2/blob/master/code/mi_functions.R
cpue <- cpue %>% 
  mutate(T_K = temp + 273.15,
         oxy2 = oxy * (10^3) / 22.391,
         size_cl = ifelse(length_cm > 30, "large", "small"),
         phi = ifelse(size_cl == "large", 
                      A_conc*(5000^n)*(oxy2 / exp(-E_conc / (kb * T_K))),
                      A_conc*(30^n)*(oxy2 / exp(-E_conc / (kb * T_K)))))
#> mutate: new variable 'T_K' (double) with 8,652 unique values and <1% NA
#>         new variable 'oxy2' (double) with 8,651 unique values and <1% NA
#>         new variable 'size_cl' (character) with 2 unique values and 0% NA
#>         new variable 'phi' (double) with 17,303 unique values and <1% NA

ggplot(cpue, aes(phi)) + geom_histogram() + facet_wrap(~size_cl, ncol = 1)
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 4487 rows containing non-finite values (stat_bin).


cpue <- cpue %>%
  mutate(
         # log_depth = log(depth),
         # log_depth2 = log_depth*log_depth,
         # depth_sc = (log_depth - mean_depth) / sd_depth,
         # depth_sc_sq = (log_depth2 - mean_depth2) / sd_depth2,
         depth_sc = (depth - mean_depth) / sd_depth,
         temp_sc = (temp - mean_temp) / sd_temp,
         oxy_sc = (oxy - mean_oxy) / sd_oxy,
         phi_sc = (phi - mean_phi) / sd_phi)
#> mutate: new variable 'depth_sc' (double) with 150 unique values and 0% NA
#>         new variable 'temp_sc' (double) with 8,652 unique values and <1% NA
#>         new variable 'oxy_sc' (double) with 8,651 unique values and <1% NA
#>         new variable 'phi_sc' (double) with 17,303 unique values and <1% NA

Summarise data

# Summaries density by length group per haul
cod_cpue <- cpue %>%
  filter(species == "cod") %>%
  mutate(size_group = ifelse(length_cm >= 30, "large", "small")) %>% 
  group_by(haul_id, size_group) %>% 
  summarise(density = sum(density)) %>% 
  ungroup()
#> filter: removed 539,125 rows (32%), 1,156,224 rows remaining
#> mutate: new variable 'size_group' (character) with 2 unique values and 0% NA
#> group_by: 2 grouping variables (haul_id, size_group)
#> summarise: now 18,066 rows and 3 columns, one group variable remaining (haul_id)
#> ungroup: no grouping variables

  cod_cpue_small <- cod_cpue %>%
    filter(size_group == "small") %>%
    dplyr::select(-size_group) %>% 
    rename(density_cod_small = density)
#> filter: removed 9,033 rows (50%), 9,033 rows remaining
#> rename: renamed one variable (density_cod_small)

  cod_cpue_large <- cod_cpue %>%
    filter(size_group == "large") %>%
    dplyr::select(-size_group) %>% 
    rename(density_cod_large = density)
#> filter: removed 9,033 rows (50%), 9,033 rows remaining
#> rename: renamed one variable (density_cod_large)

fle_cpue <- cpue %>%
  filter(species == "flounder") %>%
  mutate(size_group = ifelse(length_cm >= 22, "large", "small")) %>% 
  group_by(haul_id, size_group) %>% 
  summarise(density = sum(density)) %>% 
  ungroup()
#> filter: removed 1,156,224 rows (68%), 539,125 rows remaining
#> mutate: new variable 'size_group' (character) with 2 unique values and 0% NA
#> group_by: 2 grouping variables (haul_id, size_group)
#> summarise: now 17,622 rows and 3 columns, one group variable remaining (haul_id)
#> ungroup: no grouping variables

  fle_cpue_small <- fle_cpue %>%
    filter(size_group == "small") %>%
    dplyr::select(-size_group) %>% 
    rename(density_fle_small = density)
#> filter: removed 8,811 rows (50%), 8,811 rows remaining
#> rename: renamed one variable (density_fle_small)

  fle_cpue_large <- fle_cpue %>%
    filter(size_group == "large") %>%
    dplyr::select(-size_group) %>% 
    rename(density_fle_large = density)
#> filter: removed 8,811 rows (50%), 8,811 rows remaining
#> rename: renamed one variable (density_fle_large)


# Check 0 catches
cod_cpue_small %>%
  group_by(haul_id) %>%
  summarise(haul_dens = sum(density_cod_small)) %>%
  ungroup() %>%
  filter(!haul_dens == 0)
#> group_by: one grouping variable (haul_id)
#> summarise: now 9,033 rows and 2 columns, ungrouped
#> ungroup: no grouping variables
#> filter: removed 1,685 rows (19%), 7,348 rows remaining
#> # A tibble: 7,348 × 2
#>    haul_id                  haul_dens
#>    <chr>                        <dbl>
#>  1 1993:1:GFR:SOL:H20:21:1   15.1    
#>  2 1993:1:GFR:SOL:H20:23:31   0.00919
#>  3 1993:1:GFR:SOL:H20:25:2  278.     
#>  4 1993:1:GFR:SOL:H20:26:3  112.     
#>  5 1993:1:GFR:SOL:H20:27:27   3.64   
#>  6 1993:1:GFR:SOL:H20:28:24 741.     
#>  7 1993:1:GFR:SOL:H20:29:29   0.174  
#>  8 1993:1:GFR:SOL:H20:31:25  20.3    
#>  9 1993:1:GFR:SOL:H20:32:5  571.     
#> 10 1993:1:GFR:SOL:H20:33:6  282.     
#> # … with 7,338 more rows

cod_cpue_large %>%
  group_by(haul_id) %>%
  summarise(haul_dens = sum(density_cod_large)) %>%
  ungroup() %>%
  filter(!haul_dens == 0)
#> group_by: one grouping variable (haul_id)
#> summarise: now 9,033 rows and 2 columns, ungrouped
#> ungroup: no grouping variables
#> filter: removed 1,200 rows (13%), 7,833 rows remaining
#> # A tibble: 7,833 × 2
#>    haul_id                  haul_dens
#>    <chr>                        <dbl>
#>  1 1993:1:GFR:SOL:H20:21:1     142.  
#>  2 1993:1:GFR:SOL:H20:22:32     10.3 
#>  3 1993:1:GFR:SOL:H20:23:31      5.22
#>  4 1993:1:GFR:SOL:H20:24:30     19.4 
#>  5 1993:1:GFR:SOL:H20:25:2     531.  
#>  6 1993:1:GFR:SOL:H20:26:3     251.  
#>  7 1993:1:GFR:SOL:H20:27:27     19.5 
#>  8 1993:1:GFR:SOL:H20:28:24   1561.  
#>  9 1993:1:GFR:SOL:H20:29:29     30.6 
#> 10 1993:1:GFR:SOL:H20:30:28      6.84
#> # … with 7,823 more rows

fle_cpue_small %>%
  group_by(haul_id) %>%
  summarise(haul_dens = sum(density_fle_small)) %>%
  ungroup() %>%
  filter(!haul_dens == 0)
#> group_by: one grouping variable (haul_id)
#> summarise: now 8,811 rows and 2 columns, ungrouped
#> ungroup: no grouping variables
#> filter: removed 3,236 rows (37%), 5,575 rows remaining
#> # A tibble: 5,575 × 2
#>    haul_id                  haul_dens
#>    <chr>                        <dbl>
#>  1 1993:1:GFR:SOL:H20:21:1      0.554
#>  2 1993:1:GFR:SOL:H20:23:31     0.382
#>  3 1993:1:GFR:SOL:H20:26:3      1.16 
#>  4 1993:1:GFR:SOL:H20:27:27     3.50 
#>  5 1993:1:GFR:SOL:H20:29:29     0.498
#>  6 1993:1:GFR:SOL:H20:30:28     0.297
#>  7 1993:1:GFR:SOL:H20:31:25     0.526
#>  8 1993:1:GFR:SOL:H20:32:5      0.779
#>  9 1993:1:GFR:SOL:H20:33:6      1.51 
#> 10 1993:1:GFR:SOL:H20:34:7     10.1  
#> # … with 5,565 more rows

fle_cpue_large %>%
  group_by(haul_id) %>%
  summarise(haul_dens = sum(density_fle_large)) %>%
  ungroup() %>%
  filter(!haul_dens == 0)
#> group_by: one grouping variable (haul_id)
#> summarise: now 8,811 rows and 2 columns, ungrouped
#> ungroup: no grouping variables
#> filter: removed 1,024 rows (12%), 7,787 rows remaining
#> # A tibble: 7,787 × 2
#>    haul_id                  haul_dens
#>    <chr>                        <dbl>
#>  1 1993:1:GFR:SOL:H20:21:1      18.3 
#>  2 1993:1:GFR:SOL:H20:22:32     13.3 
#>  3 1993:1:GFR:SOL:H20:23:31      1.52
#>  4 1993:1:GFR:SOL:H20:24:30      3.78
#>  5 1993:1:GFR:SOL:H20:25:2      33.4 
#>  6 1993:1:GFR:SOL:H20:26:3      25.8 
#>  7 1993:1:GFR:SOL:H20:27:27     20.9 
#>  8 1993:1:GFR:SOL:H20:28:24     21.2 
#>  9 1993:1:GFR:SOL:H20:29:29      7.88
#> 10 1993:1:GFR:SOL:H20:30:28      9.28
#> # … with 7,777 more rows
cod_cpue_small <- left_join(cod_cpue_small, 
                            cpue %>%
                              filter(species == "cod") %>% 
                              distinct(haul_id, .keep_all = TRUE) %>%
                              dplyr::select(haul_id, year, quarter, X, Y, oxy, temp, depth,
                                             oxy_sc, temp_sc, depth_sc, phi, phi_sc)) %>% 
  drop_na(phi)
#> filter: removed 539,125 rows (32%), 1,156,224 rows remaining
#> distinct: removed 1,147,191 rows (99%), 9,033 rows remaining
#> Joining, by = "haul_id"
#> left_join: added 12 columns (year, quarter, X, Y, oxy, …)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 9,033
#> > =======
#> > rows total 9,033
#> drop_na: removed 26 rows (<1%), 9,007 rows remaining

cod_cpue_large <- left_join(cod_cpue_large,
                            cpue %>%
                              filter(species == "cod") %>% 
                              distinct(haul_id, .keep_all = TRUE) %>%
                              dplyr::select(haul_id, year, quarter, X, Y, oxy, temp, depth,
                                             oxy_sc, temp_sc, depth_sc, phi, phi_sc)) %>% 
  drop_na(phi)
#> filter: removed 539,125 rows (32%), 1,156,224 rows remaining
#> distinct: removed 1,147,191 rows (99%), 9,033 rows remaining
#> Joining, by = "haul_id"left_join: added 12 columns (year, quarter, X, Y, oxy, …)
#>            > rows only in x       0
#>            > rows only in y  (    0)
#>            > matched rows     9,033
#>            >                 =======
#>            > rows total       9,033
#> drop_na: removed 26 rows (<1%), 9,007 rows remaining

fle_cpue_small <- left_join(fle_cpue_small, 
                            cpue %>%
                              filter(species == "flounder") %>% 
                              distinct(haul_id, .keep_all = TRUE) %>%
                              dplyr::select(haul_id, year, quarter, X, Y, oxy, temp, depth,
                                             oxy_sc, temp_sc, depth_sc))
#> filter: removed 1,156,224 rows (68%), 539,125 rows remaining
#> distinct: removed 530,314 rows (98%), 8,811 rows remaining
#> Joining, by = "haul_id"left_join: added 10 columns (year, quarter, X, Y, oxy, …)
#>            > rows only in x       0
#>            > rows only in y  (    0)
#>            > matched rows     8,811
#>            >                 =======
#>            > rows total       8,811

fle_cpue_large <- left_join(fle_cpue_large,
                            cpue %>%
                              filter(species == "flounder") %>% 
                              distinct(haul_id, .keep_all = TRUE) %>%
                              dplyr::select(haul_id, year, quarter, X, Y, oxy, temp, depth,
                                            oxy_sc, temp_sc, depth_sc))
#> filter: removed 1,156,224 rows (68%), 539,125 rows remaining
#> distinct: removed 530,314 rows (98%), 8,811 rows remaining
#> Joining, by = "haul_id"left_join: added 10 columns (year, quarter, X, Y, oxy, …)
#>            > rows only in x       0
#>            > rows only in y  (    0)
#>            > matched rows     8,811
#>            >                 =======
#>            > rows total       8,811

Fit models

# Small cod
spde_small_cod <- make_mesh(cod_cpue_small, xy_cols = c("X", "Y"),
                            n_knots = 150, 
                            type = "kmeans", seed = 42)


# Large cod
spde_large_cod <- make_mesh(cod_cpue_large, xy_cols = c("X", "Y"),
                            n_knots = 150, 
                            type = "kmeans", seed = 42)

# Small flounder
spde_small_fle <- make_mesh(fle_cpue_small, xy_cols = c("X", "Y"),
                            n_knots = 150, 
                            type = "kmeans", seed = 42)


# Large flounder
spde_large_fle <- make_mesh(fle_cpue_large, xy_cols = c("X", "Y"),
                            n_knots = 150, 
                            type = "kmeans", seed = 42)

plot(spde_small_cod)
plot(spde_large_cod)

plot(spde_small_fle)
plot(spde_large_fle)

Small cod

# Small cod model
mcod_s <- sdmTMB(density_cod_small ~ 0 + as.factor(quarter) + as.factor(year) + s(depth_sc) + oxy_sc + temp_sc,
                 data = cod_cpue_small,
                 mesh = spde_small_cod,
                 family = tweedie(link = "log"),
                 spatiotemporal = "off",
                 spatial = "on",
                 time = "year",
                 reml = FALSE)
#> Warning: The model may not have converged. Maximum final gradient:
#> 0.056067183534736.

mcod_s_v2 <- run_extra_optimization(mcod_s, nlminb_loops = 1, newton_loops = 1)
sanity(mcod_s_v2)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No fixed-effect standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
# Small cod model
mcod_s1 <- sdmTMB(density_cod_small ~ 0 + as.factor(quarter) + as.factor(year) + s(depth_sc) + breakpt(oxy_sc) + temp_sc,
                 data = cod_cpue_small,
                 mesh = spde_small_cod,
                 family = tweedie(link = "log"),
                 spatiotemporal = "off",
                 spatial = "on",
                 time = "year",
                 reml = FALSE)
#> Warning: The model may not have converged. Maximum final gradient:
#> 0.011874176036468.

mcod_s1_v2 <- run_extra_optimization(mcod_s1, nlminb_loops = 1, newton_loops = 1)
sanity(mcod_s1_v2)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No fixed-effect standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
# Small cod model
mcod_s2 <- sdmTMB(density_cod_small ~ 0 + as.factor(quarter) + as.factor(year) + s(depth_sc) + oxy_sc * temp_sc,
                 data = cod_cpue_small,
                 mesh = spde_small_cod,
                 family = tweedie(link = "log"),
                 spatiotemporal = "off",
                 spatial = "on",
                 time = "year",
                 reml = FALSE)
#> Warning: The model may not have converged. Maximum final gradient:
#> 0.0247736933054483.

mcod_s2_v2 <- run_extra_optimization(mcod_s2, nlminb_loops = 1, newton_loops = 1)
sanity(mcod_s2_v2)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
# Small cod model
mcod_s3 <- sdmTMB(density_cod_small ~ 0 + as.factor(quarter) + as.factor(year) + s(depth_sc) + breakpt(phi_sc),
                 data = cod_cpue_small,
                 mesh = spde_small_cod,
                 family = tweedie(link = "log"),
                 spatiotemporal = "off",
                 spatial = "on",
                 time = "year",
                 reml = FALSE)

mcod_s3_v2 <- run_extra_optimization(mcod_s3, nlminb_loops = 1, newton_loops = 1)
sanity(mcod_s3_v2)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No fixed-effect standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
# Small cod model
mcod_s4 <- sdmTMB(density_cod_small ~ 0 + as.factor(quarter) + as.factor(year) + s(depth_sc) + s(phi_sc),
                 data = cod_cpue_small,
                 mesh = spde_small_cod,
                 family = tweedie(link = "log"),
                 spatiotemporal = "off",
                 spatial = "on",
                 time = "year",
                 reml = FALSE)
#> Warning: The model may not have converged. Maximum final gradient:
#> 0.0109864218642901.

mcod_s4_v2 <- run_extra_optimization(mcod_s4, nlminb_loops = 1, newton_loops = 1)
sanity(mcod_s4_v2)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No fixed-effect standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#mcmc_res_sc <- residuals(mcod_s_v2, type = "mle-mcmc", mcmc_iter = 201, mcmc_warmup = 200)
#qqnorm(mcmc_res_sc);qqline(mcmc_res_sc)

Large cod

# Large cod model
mcod_l <- sdmTMB(density_cod_large ~ 0 + as.factor(quarter) + as.factor(year) + s(depth_sc) + oxy_sc + temp_sc,
                 data = cod_cpue_large,
                 mesh = spde_large_cod,
                 family = tweedie(link = "log"),
                 spatiotemporal = "off",
                 spatial = "on",
                 time = "year",
                 reml = FALSE)
#> Warning: The model may not have converged. Maximum final gradient:
#> 0.0333087100107043.

mcod_l_v2 <- run_extra_optimization(mcod_l, nlminb_loops = 1, newton_loops = 1)
sanity(mcod_l_v2)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
# Large cod model
mcod_l1 <- sdmTMB(density_cod_large ~ 0 + as.factor(quarter) + as.factor(year) + s(depth_sc) + breakpt(oxy_sc) + temp_sc,
                 data = cod_cpue_large,
                 mesh = spde_large_cod,
                 family = tweedie(link = "log"),
                 spatiotemporal = "off",
                 spatial = "on",
                 time = "year",
                 reml = FALSE)
#> Warning: The model may not have converged. Maximum final gradient:
#> 0.0103058445924128.

mcod_l1_v2 <- run_extra_optimization(mcod_l1, nlminb_loops = 1, newton_loops = 1)
sanity(mcod_l1_v2)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
# Large cod model
mcod_l2 <- sdmTMB(density_cod_large ~ 0 + as.factor(quarter) + as.factor(year) + s(depth_sc) + oxy_sc*temp_sc,
                 data = cod_cpue_large,
                 mesh = spde_large_cod,
                 family = tweedie(link = "log"),
                 spatiotemporal = "off",
                 spatial = "on",
                 time = "year",
                 reml = FALSE)
#> Warning: The model may not have converged. Maximum final gradient:
#> 0.0611529616991593.

mcod_l2_v2 <- run_extra_optimization(mcod_l2, nlminb_loops = 1, newton_loops = 1)
sanity(mcod_l2_v2)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No fixed-effect standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
# Large cod model
mcod_l3 <- sdmTMB(density_cod_large ~ 0 + as.factor(quarter) + as.factor(year) + s(depth_sc) + breakpt(phi_sc),
                 data = cod_cpue_large,
                 mesh = spde_large_cod,
                 family = tweedie(link = "log"),
                 spatiotemporal = "off",
                 spatial = "on",
                 time = "year",
                 reml = FALSE)
#> Warning: The model may not have converged. Maximum final gradient:
#> 0.0155203509221983.

mcod_l3_v2 <- run_extra_optimization(mcod_l3, nlminb_loops = 1, newton_loops = 1)
sanity(mcod_l3_v2)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No fixed-effect standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
# Large cod model
mcod_l4 <- sdmTMB(density_cod_large ~ 0 + as.factor(quarter) + as.factor(year) + s(depth_sc) + s(phi_sc),
                 data = cod_cpue_large,
                 mesh = spde_large_cod,
                 family = tweedie(link = "log"),
                 spatiotemporal = "off",
                 spatial = "on",
                 time = "year",
                 reml = FALSE)

mcod_l4_v2 <- run_extra_optimization(mcod_l4, nlminb_loops = 1, newton_loops = 1)
sanity(mcod_l4_v2)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No fixed-effect standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#mcmc_res_lc <- residuals(mcod_l_v2, type = "mle-mcmc", mcmc_iter = 201, mcmc_warmup = 200)
#qqnorm(mcmc_res_lc);qqline(mcmc_res_lc)

Small flounder

# Small flounder model
# mfle_s <- sdmTMB(density_fle_small ~ 0 + as.factor(quarter) + as.factor(year) + s(depth_sc) +
#                    breakpt(oxy_sc) + temp_sc,
#                  data = fle_cpue_small,
#                  mesh = spde_small_fle,
#                  family = tweedie(link = "log"),
#                  spatiotemporal = "off",
#                  spatial = "on",
#                  time = "year",
#                  reml = FALSE)
# 
# mfle_s_v2 <- run_extra_optimization(mfle_s, nlminb_loops = 1, newton_loops = 1)
# Small flounder model
# mfle_s2 <- sdmTMB(density_fle_small ~ 0 + as.factor(quarter) + as.factor(year) + s(depth_sc) +
#                    breakpt(oxy_sc) + temp_sc,
#                  data = fle_cpue_small,
#                  mesh = spde_small_fle,
#                  family = tweedie(link = "log"),
#                  spatiotemporal = "off",
#                  spatial = "on",
#                  time = "year",
#                  reml = FALSE)
# 
# mfle_s2_v2 <- run_extra_optimization(mfle_s2, nlminb_loops = 1, newton_loops = 1)
#mcmc_res_sf <- residuals(mfle_s_v2, type = "mle-mcmc", mcmc_iter = 201, mcmc_warmup = 200)
#qqnorm(mcmc_res_sf);qqline(mcmc_res_sf)

Large flounder

# Large flounder model
# mfle_l <- sdmTMB(density_fle_large ~ 0 + as.factor(quarter) + as.factor(year) + s(depth_sc) +
#                    breakpt(oxy_sc) + temp_sc,
#                  data = fle_cpue_large,
#                  mesh = spde_large_fle,
#                  family = tweedie(link = "log"),
#                  spatiotemporal = "off",
#                  spatial = "on",
#                  time = "year",
#                  reml = FALSE)
# 
# mfle_l_v2 <- run_extra_optimization(mfle_l, nlminb_loops = 1, newton_loops = 1)
# Large flounder model
# mfle_l2 <- sdmTMB(density_fle_large ~ 0 + as.factor(quarter) + as.factor(year) + s(depth_sc) +
#                    oxy_sc*temp_sc,
#                  data = fle_cpue_large,
#                  mesh = spde_large_fle,
#                  family = tweedie(link = "log"),
#                  spatiotemporal = "off",
#                  spatial = "on",
#                  time = "year",
#                  reml = FALSE)
# 
# mfle_l2_v2 <- run_extra_optimization(mfle_l2, nlminb_loops = 1, newton_loops = 1)
#mcmc_res_lf <- residuals(mfle_l_v2, type = "mle-mcmc", mcmc_iter = 201, mcmc_warmup = 200)
#qqnorm(mcmc_res_lf);qqline(mcmc_res_lf)

Results

AIC

AIC(mcod_s_v2, mcod_s1_v2, mcod_s2_v2, mcod_s3_v2, mcod_s4_v2)
#>            df      AIC
#> mcod_s_v2  36 81998.42
#> mcod_s1_v2 37 81908.64
#> mcod_s2_v2 37 81832.75
#> mcod_s3_v2 36 82128.92
#> mcod_s4_v2 36 81727.10
AIC(mcod_l_v2, mcod_l1_v2, mcod_l2_v2, mcod_l3_v2, mcod_l4_v2)
#>            df      AIC
#> mcod_l_v2  36 110363.0
#> mcod_l1_v2 37 110283.7
#> mcod_l2_v2 37 110213.7
#> mcod_l3_v2 36 110344.2
#> mcod_l4_v2 36 110023.9

For both small and large cod, the smooth phi model is best

ggplot(cod_cpue_large, aes(phi, density_cod_large)) + geom_point() + stat_smooth()
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

ggplot(cod_cpue_small, aes(phi, density_cod_small)) + geom_point() + stat_smooth()
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'


tidy(mcod_l4_v2)
#>                   term    estimate std.error
#> 1  as.factor(quarter)1  4.42736002 0.4252533
#> 2  as.factor(quarter)4  4.44623473 0.4251376
#> 3  as.factor(year)1994  0.17447088 0.1251899
#> 4  as.factor(year)1995  0.36223149 0.1256474
#> 5  as.factor(year)1996  0.49711011 0.1235093
#> 6  as.factor(year)1997 -0.52187125 0.1276621
#> 7  as.factor(year)1998 -0.64041354 0.1249065
#> 8  as.factor(year)1999 -0.57771772 0.1255800
#> 9  as.factor(year)2000 -0.33946949 0.1319292
#> 10 as.factor(year)2001 -0.25852717 0.1196450
#> 11 as.factor(year)2002  0.46512739 0.1191739
#> 12 as.factor(year)2003 -0.10063961 0.1235872
#> 13 as.factor(year)2004 -0.42438061 0.1231739
#> 14 as.factor(year)2005  0.22425024 0.1165817
#> 15 as.factor(year)2006  0.32941648 0.1173735
#> 16 as.factor(year)2007  0.33068053 0.1143035
#> 17 as.factor(year)2008  0.65553792 0.1138592
#> 18 as.factor(year)2009  0.77452186 0.1116785
#> 19 as.factor(year)2010  1.20435709 0.1121921
#> 20 as.factor(year)2011  0.76709423 0.1130294
#> 21 as.factor(year)2012  0.51482769 0.1158932
#> 22 as.factor(year)2013  0.21750426 0.1141133
#> 23 as.factor(year)2014  0.15737085 0.1158155
#> 24 as.factor(year)2015  0.35017962 0.1137032
#> 25 as.factor(year)2016  0.18424103 0.1133089
#> 26 as.factor(year)2017 -0.05126881 0.1139330
#> 27 as.factor(year)2018 -0.37601593 0.1156271
#> 28 as.factor(year)2019 -0.72588172 0.1387845
tidy(mcod_s4_v2)
#>                   term    estimate std.error
#> 1  as.factor(quarter)1  2.85042323 0.3406691
#> 2  as.factor(quarter)4  2.63170504 0.3412020
#> 3  as.factor(year)1994 -0.69689551 0.1593284
#> 4  as.factor(year)1995  0.12906713 0.1539344
#> 5  as.factor(year)1996 -0.29763762 0.1581533
#> 6  as.factor(year)1997 -0.86644518 0.1553678
#> 7  as.factor(year)1998 -0.44999953 0.1470866
#> 8  as.factor(year)1999 -0.45410484 0.1492065
#> 9  as.factor(year)2000 -0.13818808 0.1546355
#> 10 as.factor(year)2001  0.26872266 0.1408705
#> 11 as.factor(year)2002  0.66548736 0.1423063
#> 12 as.factor(year)2003 -0.49700848 0.1527587
#> 13 as.factor(year)2004  0.29538027 0.1437326
#> 14 as.factor(year)2005  0.69840669 0.1375207
#> 15 as.factor(year)2006  0.05817676 0.1439762
#> 16 as.factor(year)2007  0.37674494 0.1384376
#> 17 as.factor(year)2008  0.54291810 0.1375249
#> 18 as.factor(year)2009  0.54721732 0.1363738
#> 19 as.factor(year)2010  0.79860508 0.1388352
#> 20 as.factor(year)2011  0.36139637 0.1395769
#> 21 as.factor(year)2012  0.28997785 0.1413632
#> 22 as.factor(year)2013  1.04934351 0.1340687
#> 23 as.factor(year)2014  0.77760213 0.1368409
#> 24 as.factor(year)2015  0.43911386 0.1377468
#> 25 as.factor(year)2016  0.24854388 0.1372006
#> 26 as.factor(year)2017  0.27336163 0.1360151
#> 27 as.factor(year)2018  0.02334050 0.1381793
#> 28 as.factor(year)2019  0.22814807 0.1583841

large_pars <- tidy(mcod_l2_v2,
                effects = "fixed", conf.int = TRUE) %>%
  filter(!grepl('year', term)) %>%
  filter(!grepl('quarter', term)) %>% 
  mutate(size_cl = "Adult")
#> filter: removed 26 rows (84%), 5 rows remaining
#> filter: removed 2 rows (40%), 3 rows remaining
#> mutate: new variable 'size_cl' (character) with one unique value and 0% NA

small_pars <- tidy(mcod_s2_v2,
                effects = "fixed", conf.int = TRUE) %>%
  filter(!grepl('year', term)) %>%
  filter(!grepl('quarter', term)) %>% 
  mutate(size_cl = "Juvenile")
#> filter: removed 26 rows (84%), 5 rows remaining
#> filter: removed 2 rows (40%), 3 rows remaining
#> mutate: new variable 'size_cl' (character) with one unique value and 0% NA

coefs <- bind_rows(large_pars, small_pars)

coefs <- coefs %>% 
    mutate(term2 = term,
           term2 = ifelse(term == "temp_sc", "Temperature", term2),
           term2 = ifelse(term == "oxy_sc", "Oxygen", term2),
           term2 = ifelse(term == "depth_sc", "Depth", term2),
           term2 = ifelse(term == "oxy_sc:temp_sc", "Temperature x Oxygen", term2)
           )
#> mutate: new variable 'term2' (character) with 3 unique values and 0% NA

coefs <- coefs %>% 
    mutate(term2 = term,
           term2 = ifelse(term == "temp_sc", "Temperature", term2),
           term2 = ifelse(term == "oxy_sc", "Oxygen", term2),
           term2 = ifelse(term == "depth_sc", "Depth", term2),
           term2 = ifelse(term == "oxy_sc:temp_sc", "Temperature x Oxygen", term2)
           )
#> mutate: no changes

# Coefficients
ggplot(coefs, aes(term2, estimate, color = size_cl)) +
  geom_hline(yintercept = 0, linetype = 2, color = "gray50", size = 0.5) +
  geom_point(size = 3, position = position_dodge(width = 0.2)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0, size = 1,
                position = position_dodge(width = 0.2)) +
  labs(x = "Predictor", y = "Standardized coefficient") +
  theme(legend.position = "right",
        axis.text.x = element_text(angle = 90)) +
  NULL 

2. Marginal effects plot

Weird that it’s unimodal and declines at high phi’s. Maybe super coastal?

# Marginal
nd <- data.frame(
  phi_sc = seq(min(cod_cpue_small$phi_sc) + 0.5,
               max(cod_cpue_small$phi_sc) - 0.2, length.out = 100), 
  year = 2015L,
  quarter = 1,
  oxy_sc = mean(cod_cpue_small$oxy_sc),
  temp_sc = mean(cod_cpue_small$temp_sc),
  depth_sc = mean(cod_cpue_small$depth_sc)
)

p_adu <- predict(mcod_l4_v2, newdata = nd, se_fit = TRUE, re_form = NA, xy_cols = c("X", "Y"))
p_juv <- predict(mcod_s4_v2, newdata = nd, se_fit = TRUE, re_form = NA, xy_cols = c("X", "Y"))

p <- bind_rows(p_adu %>% mutate(size_cl = "Adult"),
               p_juv %>% mutate(size_cl = "Juvenile")) %>% 
  mutate(phi = (phi_sc * sd_phi) + mean_phi)
#> mutate: new variable 'size_cl' (character) with one unique value and 0% NA
#> mutate: new variable 'size_cl' (character) with one unique value and 0% NA
#> mutate: new variable 'phi' (double) with 100 unique values and 0% NA

ggplot(p, aes(phi, exp(est), color = size_cl, fill = size_cl, 
              ymin = exp(est - 1.96 * est_se), 
              ymax = exp(est + 1.96 * est_se))) +
  geom_line() +
  geom_ribbon(alpha = 0.4, color = NA) + 
  geom_vline(xintercept = 1, linetype = 2, color = "grey30") +
  coord_cartesian(expand = 0) + 
  scale_x_continuous(breaks = seq(-2, 40, by = 2))

  1. Predictions on maps

pred_grid
#> # A tibble: 249,453 × 14
#>        X     Y depth  year   oxy  temp   lon   lat sub_div ices_…¹ depth…²   T_K
#>    <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>   <dbl> <chr>     <dbl> <dbl>
#>  1   306  6012  13.2  1993  4.11  9.78  12.0  54.2      24 37G2         21  283.
#>  2   306  6012  13.2  1994  6.27  9.66  12.0  54.2      24 37G2         21  283.
#>  3   306  6012  13.2  1995  5.42  9.63  12.0  54.2      24 37G2         21  283.
#>  4   306  6012  13.2  1996  6.11  9.99  12.0  54.2      24 37G2         21  283.
#>  5   306  6012  13.2  1997  5.78  9.84  12.0  54.2      24 37G2         21  283.
#>  6   306  6012  13.2  1998  6.45  8.56  12.0  54.2      24 37G2         21  282.
#>  7   306  6012  13.2  1999  5.53 10.2   12.0  54.2      24 37G2         21  283.
#>  8   306  6012  13.2  2000  4.37 12.5   12.0  54.2      24 37G2         21  286.
#>  9   306  6012  13.2  2001  5.81  9.91  12.0  54.2      24 37G2         21  283.
#> 10   306  6012  13.2  2002  5.26  9.92  12.0  54.2      24 37G2         21  283.
#> # … with 249,443 more rows, 2 more variables: oxy2 <dbl>, phi <dbl>, and
#> #   abbreviated variable names ¹​ices_rect, ²​depth_sd

plot_map_fc + 
  geom_raster(data = pred_grid, aes(X*1000, Y*1000, fill = phi)) + 
  facet_wrap(~year) + 
  scale_fill_viridis()
#> Warning: Removed 9882 rows containing missing values (geom_raster).

  1. weighted phi by size group and year
pred_grid <- pred_grid %>%
  mutate(quarter = factor(4),
         depth_sc = (depth - mean(depth)) / sd(depth),
         phi_sc = (phi - mean(phi)) / sd(phi))
#> mutate: new variable 'quarter' (factor) with one unique value and 0% NA
#>         new variable 'depth_sc' (double) with 6,238 unique values and 0% NA
#>         new variable 'phi_sc' (double) with 249,453 unique values and 0% NA

large_pred <- predict(mcod_l4_v2, newdata = pred_grid)
small_pred <- predict(mcod_s4_v2, newdata = pred_grid)

# Small cod
wm_phi_small <- small_pred %>%
  group_by(year) %>%
  summarise("Median" = hutils::weighted_quantile(v = phi, w = exp(est), p = c(0.5)),
            "1st quartile" = hutils::weighted_quantile(v = phi, w = exp(est), p = c(0.25)),
            "3rd quartile" = hutils::weighted_quantile(v = phi, w = exp(est), p = c(0.75))) %>% 
  ungroup() %>% 
  mutate(group = "Juvenile") 
#> group_by: one grouping variable (year)
#> summarise: now 27 rows and 4 columns, ungrouped
#> ungroup: no grouping variables
#> mutate: new variable 'group' (character) with one unique value and 0% NA

# Large cod
wm_phi_large <- large_pred %>%
  group_by(year) %>%
  summarise("Median" = hutils::weighted_quantile(v = phi, w = exp(est), p = c(0.5)),
            "1st quartile" = hutils::weighted_quantile(v = phi, w = exp(est), p = c(0.25)),
            "3rd quartile" = hutils::weighted_quantile(v = phi, w = exp(est), p = c(0.75))) %>% 
  ungroup() %>% 
  mutate(group = "Adult")
#> group_by: one grouping variable (year)
#> summarise: now 27 rows and 4 columns, ungrouped
#> ungroup: no grouping variables
#> mutate: new variable 'group' (character) with one unique value and 0% NA

env_phi <- pred_grid %>%
  group_by(year) %>%
  summarise("Median" = quantile(phi, p = c(0.5)),
            "1st quartile" = quantile(phi, p = c(0.25)),
            "3rd quartile" = quantile(phi, p = c(0.75))) %>% 
  ungroup() %>% 
  mutate(group = "Environment") 
#> group_by: one grouping variable (year)
#> summarise: now 27 rows and 4 columns, ungrouped
#> ungroup: no grouping variables
#> mutate: new variable 'group' (character) with one unique value and 0% NA
  
phi_series <- bind_rows(wm_phi_small, wm_phi_large, env_phi)

ggplot(phi_series, aes(year, Median, color = group, fill = group)) +
  geom_line() +
  geom_ribbon(aes(ymin = `1st quartile`, ymax = `3rd quartile`), alpha = 0.2, color = NA) +
  facet_wrap(~group) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 6)) +
  guides(fill = "none", color = "none") +
  labs(y = "phi", x = "Year") +
  theme_plot() +
  theme(axis.text.x = element_text(angle = 90)) +
  NULL


ggplot(phi_series, aes(year, Median, color = group)) +
  geom_line() +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 6)) +
  labs(y = "phi", x = "Year", color = "") +
  theme_plot() +
  theme(axis.text.x = element_text(angle = 90)) +
  NULL